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Summary 

A  3-D  electromagnetic  particle-in-cell  code  (EMPIC)  was  developed  for  large-scale  plasma 
simulation  on  parallel  computers.  To  simulate  plasma  problems  with  complex  geometries  such  as  high- 
power  microwave  generation  devices,  curvilinear  coordinates  were  used.  A  logically  connected  Cartesian 
grid  consisted  of  hexahedral  cells  that  could  be  deformed  to  body-fit  complex  shapes.  A  finite- volume 
method  for  a  non-orthogonal  grid  was  used  to  calculate  the  electromagnetic  field.  This  method  is  based  on 
Gedney-Lansing  [1]  and  Madsen  [2]  algorithms,  and  is  reduced  to  a  standard  FDTD  algorithm  for  an 
orthogonal  grid.  Particle  updates  use  current  deposit  formulation  of  Villasenor  and  Buneman  [3] 
generalized  to  non-orthogonal  meshes,  preserving  charge  and  current. 

The  numerical  stability  of  the  electromagnetic  algorithm  was  analyzed  for  a  planar  EM  wave 
propagation  in  a  distorted  periodic  box.  The  influence  of  face-to-edge  transformation  properties  on  the  EM 
algorithm  stability  is  analyzed.  Energy  conservation  for  a  full  PIC  simulation  of  two-stream  instability  on 
non-orthogonal  meshes  is  studied. 

Governing  equations. 


The  EM  fields  in  the  code  are  computed  from  time-dependent  Maxwell’s  equations  in  an  integral 

form: 


<9,  jE  ds  =  cjB-dl  +  Jy  ds 
d,\Bds  =  -c§E-dl 


The  particle  trajectories  are  calculated  from  relativistic  Newton’s  second  law  for  a  single  charged 
particle  in  an  electromagnetic  field: 


dt  ymV  =  q 
d,X  —  V 
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The  EM  algorithm  is  based  on  explicit  DSI  (discrete  surface  integral)  solution  of  Maxwell’s 
equations,  tailored  to  the  case  of  hexahedral  grid  cells.  Particle  trajectories  are  calculated  by  a  standard 
leap-frog  algorithm  described  in  Birdsall  and  Langdon  [4],  The  code  reduces  to  a  standard  orthogonal 
EMPIC  code  on  Cartesian  grid,  (Wang,  Liewer  and  Decyk,  [5]). 


Electromagnetic  Algorithm  Numerical  Stability  Study. 


The  algorithm  uses  a  staggered  grid  system.  The  primary  cell  B  grid  and  dual  cell  E  grid  are 
defined.  The  fundamental  code  variables  are  B-ds  (on  each  face  of  B  grid  cells),  E-ds( on  each  face  of  E 
cells),  B-dl  (on  edges  of  the  E  grid),  and  E-dl  (on  edges  of  the  B  grid),  Fig.  1 . 

The  integral  Maxwell’s  equations  are  advanced  using  a  2nd  order  leap-frog  scheme: 

/  \n+l/2 

( E  •  dsj+l  =  ( E  •  dsj  +  A t  f  cX  B  •  dl  +  J  •  ds 

V  edges  J 


/  \n+3/2  /  \n  +  l/2 

(B-ds)  =  [B  ■  ds)  -A  tc 
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J^Edl 

\  edges  J 


The  algorithm  is  completely  specified  with  transformation  from  face  quantities  (B-ds,  E-ds)  to  edge 
quantities.  For  a  uniform  grid  this  transformation  is  trivial,  so  that 


Bdl-Bds 


d l 

ds 


In  this  limit,  a  standard  2nd-order  Yee  algorithm  is  recovered.  For  a  non-orthogonal  grid,  dl  is  not 
parallel  to  ds,  so  additional  work  has  to  be  done  (Fig.  2,  for  a  2-D  case).  There,  5  B-ds  values  are  needed  to 
compute  one  B-dl  value.  To  obtain  B-dl  at  a  particular  cell  face,  say  face  1  on  Fig.  2,  one  first  needs  to  find 
four  By,  (/=  2,5)  values  at  two  vortices  of  the  face  1  by  solving  linear  system  of  equations: 

Bu  '  ds,  =  (B  ■  ds\ 

/  \  ’i=2'5- 

Bu  ■  ds,  =  (B  ■  ds\ 

where  (B-ds)i  0=1,5)  are  known  fundamental  variables,  and  Bn  are  unknown  vectors.  Next,  a  weighted 
average  of  Blt  is  dotted  with  the  dual-edge  vector  dli  to  form  B-dl: 

Here  we  can  distinguish  tree  different  methods  of  weighting  scheme: 

Method  A:  Equal  weights:  (0  .  =  0.25 


_  \dsl  x  ds{ 

Method  B:  Volume-based  weights:  (0^  = 

Method  C:  Here  B  dl  is  a  sum  along  two  line  segments  dln  and  dl,2  from  neighboring  cell  centers  to  cell 
face: 

(Bdl\  =  (Bdl)n  +  (Bdl)u 

(B  •  dl)u  —  (-^2^12  ^13^13)^11 

{B  •  dl)n  =  [Bu(Ou  +  Bl5(Ol5)dln 
dlx  =  dln  +  dln 

where  weights  are  one-sided  volume  averaged.  The  purpose  of  Method  C  is  to  compute  portions  of  edge 
values  locally  at  each  cell  so  that  no  interprocessor  communication  would  be  necessary.  However,  Method 
C  has  turned  out  to  be  the  poorest  one  in  terms  of  numerical  stability. 

Numerical  Stability  of  the  EM  Algorithm. 

The  EM  algorithm  has  been  reported  to  have  a  weak  numerical  stability  (Brandon  and  Rambo, 
|6]).  Here  the  EM  algorithm  stability  is  analyzed  on  a  test  computation  of  a  planar  EM  wave  propagation 
through  a  distorted  periodic  box.  Numerical  parameters  for  a  test  case  are  as  follows: 
the  wave  is  propagating  along  the  z-axis, 
number  of  grid  points,  Nx=  1 7,  N>=  1 7,  9; 

box  dimensions,  .£,*=10,  Ly=10,  L~ 20; 
speed  of  light  c=l, 
wave  transient  time  is  20. 

The  grid  is  2-D  distorted  in  the  x-y  plane: 

z(i,j,k)  =  k-  h, 

x(i,j,k )  =  ihx  +  d sin (mxihx)  sin (myjhy) 
y(i,j,k)  =  /  hy  +  d  sin(mjhx )  sin(mvy7zv ) 

where  h„  hy,  h:  are  mesh  sizes  for  a  uniform  orthogonal  grid, 


2k  2k 


and  d  is  distortion  parameter.  Examples  of  the  distorted  grid  for  several  values  of  d  are  shown  on  Fig.  Note 
that  in  this  geometry  the  non-zero  value  of  the  z-component  of  E-ds  represents  numerical  error,  which 
makes  it  useful  for  analytical  purposes.  Logarithmic  plots  of  total  energy  and  maximum  value  of  (E-ds) z 
over  the  computational  domain  as  a  function  of  time  are  shown  in  Fig.  4.  Three  methods  of  weighting 
scheme  A,  B  and  C  are  compared  for  a  test  run  with  d= 0.8.  All  three  methods  are  numerically  unstable; 
Method  A  has  the  smallest  growth  rate,  while  method  C  has  the  highest  one.  Method  A  takes  7500,  and 
method  C  50  wave  transient  times  to  become  unstable.  One  can  observe  that  numerical  error  in  z 
components  of  EM  field  directly  affects  numerical  stability.  Once  the  numerical  error  becomes  comparable 
with  the  non-zero  EM  field  components  of  exact  solution,  the  instability  appears  as  an  exponential  growth 
of  total  energy.  The  growth  rate  increases  with  distortion  size  (Fig.  5). 

Langdon  [7]  suggested  the  importance  of  face-to-edge  transformation  properties  for  numerical 
stability.  He  considered  a  2-D  field  in  2-D  geometry  and  found  that  for  the  numerical  stability,  the  face-to- 
edge  transformation  must  be  symmetric.  In  the  2-D  case  shown  on  Fig  2,  it  means  that  contribution  of 
(B-ds)i  to  (B  dl)i  must  be  equal  to  the  contribution  of  (B-ds),  to  (B-dl)2. 

Here  the  analysis  for  a  3-D  DSI  algorithm  is  presented.  The  discretized  Maxwell’s  equations  in  an 
integral  form  can  be  written  as 

p  =  s  ■  R 

^ DF  ^E  LJ DE 

R  =  S  -  F 

UpF  ^  PE 

where  matrices  Se  and  Sb  are  related  as 

S  =  -S  1 


and  E,)f,  BPF,  BDE,  Epe  are  vectors  of  grid  values  of  E-ds,  B  ds,  B  dl  and  E  dl,  respectively.  The  face-to- 
edge  transformation  can  be  specified  as 


BDE  -  TB  •  BPF 


'PE 


T  •  E 

lE  ^DF 


Using  these  definitions,  we  derive 


Ebf  =  -MeEdf,  Me  =  SeTbSetTe 

bpf  =  -mbbpf,  mb  =  s/tesetb 

Thus,  for  computations  to  be  stable  numerically  below  Courant  limit,  matrices  Me  and  Mb  must  be 
symmetric  and  positive  definite.  In  this  case,  all  eigenvalues  of  Me  and  Mh  are  positive.  This  requirement 
imposes  restrictions  on  both  Te  and  Tb.  In  particular,  if 

SF  T„  =  Tet  Se  .  (i) 

then  Me  and  Mb  are  symmetric.  Let  us  consider  several  special  cases. 

1 .  Uniform  Cartesian  Mesh: 

TB=TE  =  A  (diagonal  matrix) 


Clearly,  the  condition  (1)  is  satisfied  and  solution  is  stable,  which  is  known  from  the  properties  of  the  Yee 
algorithm. 

2.  Non-orthogonal  uniformly  skewed  grid  (Fig.  6): 

Here  we  have  symmetric  face-to-edge  transformation: 


Te  =  Tb  =  te 


One  can  prove  that  condition  (1)  holds  for  such  grid.  Figure  7  shows  that  the  solution  is  indeed  numerically 
stable  for  a  uniformly  skewed  grid. 

3.  General  non-orthogonal  grid,  such  as  on  Fig.  3: 

Here  condition  (1)  does  not  hold,  so  we  can  expect  the  solution  to  be  numerically  unstable.  One  can 
speculate  that  for  different  weighting  schemes  A,  B  and  C,  condition  (1)  holds  to  a  different  extent,  which 
explains  different  results  for  numerical  stability. 

4.  Langdon  2-D  Field  Analysis: 

Langdon  considered  the  case  with 

t.  =  a 


Thus,  matrix  Mb  becomes  : 


—  $e  TeSe 


so  for  its  symmetry  Te  .  must  be 


symmetric.  However,  this  is  a  special  case.  But  in  general,  condition  (1)  must  be  satisfied  for  numerical 
solution  to  be  stable.  Condition  ( 1 )  imposes  that  for  a  3-D  field,  even  a  symmetric  transformation  does  not 
guarantee  numerical  stability. 


Test  Case  for  a  Full  PIC  Simulation  on  Non-Orthogonal  Meshes. 


Let  us  consider  a  counter-streaming  beam  system:  two  equal  electron  beams  are  set  to  move  at 
different  directions  with  opposite  velocities  Vi=0.4c.  The  electrons  within  each  beam  have  Maxwellian 
distribution  with  thermal  velocity  vt*=0.05c.  The  ions  are  considered  to  be  a  fixed  background.  This 
counter-streaming  system  generates  the  well-known  two-stream  instability.  Calculation  is  again  performed 
in  a  distorted  periodic  box  as  described  earlier.  In  Fig.  8  electromagnetic  energy,  particle  kinetic  energy  and 
total  energy  in  the  system  are  plotted  as  functions  of  time  for  distorted  and  undistorted  grids.  Figure  8 
shows  that  calculations  correctly  predict  energy  transfer  between  the  particles  and  the  fields  which  results 
from  instability  excitation  and  saturation.  The  total  energy  remains  constant,  as  it  should  be.  The  error  in 
total  energy  conservation  does  not  exceed  2%  even  for  a  highly  distorted  grid  (Fig.  9),  and  it  decreases  with 
increased  number  of  particles  (Fig.  10).  The  EM  numerical  instability  plays  no  role  in  this  computation  as 
the  characteristic  time  of  the  calculation  is  only  several  transient  wave  times. 
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